Seasonal estimation of groundwater vulnerability

Index-based methods estimate a fixed value of groundwater vulnerability (GWV); however, the effects of time variations on this estimation have not been comprehensively studied. It is imperative to estimate a time-variant vulnerability that accounts for climatic changes. In this study, we used a Pesticide DRASTICL method separating hydrogeological factors into dynamic and static groups followed by correspondence analysis. The dynamic group is composed of depth and recharge, and the static group is composed of aquifer media, soil media, topography slope, impact of vadose zone, aquifer conductivity and land use. The model results were 42.25–179.89, 33.93–159.81, 34.08–168.74, and 45.56–205.20 for spring, summer, autumn, and winter, respectively. The results showed a moderate correlation between the model predictions and observed nitrogen concentrations with R2 = 0.568 and a high correlation for phosphorus concentrations with R2 = 0.706. Our results suggest that the time-variant GWV model provides a robust yet flexible method for investigating seasonal changes in GWV. This model is an improvement to the standard index-based methods, making them sensitive to climatic changes and portraying a true vulnerability estimation. Finally, the correction of the rating scale value fixes the problem of overestimation in standard models.

www.nature.com/scientificreports/ The hydrological subregions are the Hopkins River Basin (East), which includes the Hopkins and the Merry rivers; the Portland Coastal (Southcentral), which includes the Moyne, Eumeralla, Fitzroy and Surry rivers; and the Glenelg River Basin (Northcentral and West), which includes the tributaries of the Glenelg River. The largest water body in the basin is the Rocklands Reservoir 45 .
Different soil lithologies are present in the region, with the most abundant lithologies including basalts, sedimentary, duricrust, and aeolian, followed by alluvium, limestone, alluvial, granite, colluvial, lacustrine/ aeolian, and volcanic. The least abundant soils are fluvial, aeolian and lagoonal 45 . The hydrogeological features of the study area can be seen in Fig. 2. Methods. We developed a novel method comprising five steps: (1) Data preparation. (2) Calculation of the Pesticide DRASTICL vulnerability index (IDL) using dynamic factors (D and R) and static hydrogeological factors (A, S, T, I, C, and L). (3) Independent computation of the correspondence analysis (CA) for both dynamic and static factor groups. (4) Upscaling the resulting eigenvalues to the Pesticide DRASTICL weight scale (1 to 5). (5) Recalculation of the seasonal IDL using the calculated weights from step 4.
The sources of data used for Pesticide DRASTICL factors are summarized in Table 1. Both dynamic factors were estimated with models: for D, we used HydroSight 39 , and for R, we used data resulting from the groundwater transient model for the Glenelg Hopkins Catchment Management Authority (CMA). The transient model included a three-stage approach: pre-development -1985, calibration-1985 to 1994, verification-1995 to 1999 and post-development-1995.   Table 2. Weights represent how important each factor is when compared to others and range from 1 to 5, with 5 being the most significant factor and 1 the least significant factor 30 . Each factor was reclassified by assigning the corresponding rating value and converted to raster format using ArcMap10.4.1. Each factor layer was portrayed in a grid of 50 m × 50 m using a georeferenced GDA 1994 Lambert Conformal Conic projection. Ranges represent the significant media type that contributes to the pollution potential, and such values change spatially and within their own scale 30 . Every range is evaluated against each other to calculate the relative significance in relation to pollution potential 30 . In DRASTIC methods, the factor range of relative significance is evaluated on a scale of 1 (least significant) to 10 (most significant), except for the net recharge factor, which has a scale of 1 to 9.
Rating values were established using an ordinal scale; however, assigning range values from 1 to 10 to differentiate the relative importance between ranges has been found to be contrary to what can be truly comprehended, as the human brain can only categorize effects on a scale of 1 to 9, which is known as the scale of intensity 23,30,52 . Hence, standard index-based methods tend to overestimate resulting scores by using 10 different categories when in reality only 9 categories are comprehensible by the human brain. In this study, the method uses a range scale from 1 to 9 to address the aforementioned limitation. The new rating value ranges are shown in Table 2.
Depth to water table (D). Seasonal groundwater depths (D S ) are shown in Fig. 3 a, b, c and d. In this study, space-time groundwater elevation maps from Peterson and Western 47 were adopted. This approach is an extension of Costelloe et al. 53 and Peterson et al. 54 . The maps were derived at a monthly time-step from 1 January 1980 to 1 August 2014. The maps were derived using an advanced multivariate geostatistical R package named HydroMap (https:// github. com/ peter son-tim-j/ Hydro Map). The approach allows the inclusion of topographic form, land surface elevation, coastline, and the physical constraint of the land surface elevation on the water table elevation. Furthermore, the factors within these predictors and the standard kriging factors (variogram range, sill and nugget and search radius) were derived using formal maximum likelihood estimation. Input data for kriging were derived from the HydroSight (http:// peter son-tim-j. github. io/ Hydro Sight/) time-series analysis of each groundwater hydrograph 55 . Groundwater hydrograph errors and outliers were omitted and temporally interpolated to a monthly time step using Peterson and Western 47,53 . In the aquifers of the study site where the differences in depth ranges are not large enough, a new range classification is suggested in Table 3. Seasonal depths are proposed and categorized by assigning each range a rating value from 1 to 9.
Net recharge (R). This study uses recharge values from a transient groundwater model 45 using monthly recharge data. The calibrated results for the steady-state and transient model normalized RMS errors were 2.47% and 2.24%, respectively. MODFLOW inputs were provided by Biosym within the Ensym model 45 . The method used to calculate recharge is described in the Australian groundwater modelling guidelines 56 . Seasonal recharge values were plotted in ArcMap 10.4.1 and categorized using values from Table 4. The groundwater recharge factor offers significant information to the model, as water recharge is the major driving mechanism for pollution transport 57 . Additionally, it provides specific information for both spatial and temporal variability. Monthly recharge values for the period from 1987 to 2017 were considered for this analysis. Seasonal recharge maps are shown in Fig. 4.
Aquifer media (A). Aquifer media (A) serves as a conduit for aquifer pollution; the larger the grain size and the more fractures there are in the aquifer, the greater the pollution potential 23 . The aquifer medium or media are responsible for controlling the route and path length for a contaminant to flow; the larger the porosity of the aquifer is, the lower the attenuation capability, and therefore, the larger the pollution potential. Aquifer configuration materials were obtained from the Victorian Aquifer Framework (VAF) 58 . The lithological units were Net recharge R 4 1-9 1-9 Aquifer media A 3 1-10 1-9 Soil media S 5 1-10 1-9 Topography slope T 3 1-10 1-9 Impact of vadose zone I 4 1-10 1-9 Hydraulic conductivity C 2 1-10 1-9 Land Use L 5 1-10 1-9 www.nature.com/scientificreports/ categorized with new rating values from 1 to 9 and are presented in Table 5. Next, values were modified from Aller et al. 23,43 and adapted to the Victoria geological units (see Fig. 5a).

Soil media (S).
Soil media is the upper portion of the vadose zone; it has significant biological activity, and it is considered the upper weathered zone of the Earth, averaging six feet (1.8 m) or less 23 . Soil media control recharge and contaminant attenuation processes such as filtration, biodegradation, sorption, and volatilization 23 . The upper portion of the soil configuration is used for this factor. Data were extracted from the Atlas of Australian Soil (ASRI), and the variable that represents the topmost layer is (TEXT_TOP) within the Australian Soil Atlas Dataset. The categorization of the Australian soil texture is found in the glossary of soil terms and in  www.nature.com/scientificreports/ "Estimation of Soil Properties using the Atlas of Australian Soils" 59 . A rating for soil texture groups is presented in Table 6, and the map for the soil media is shown in Fig. 5c.
Topography slope (T). Topography refers to the slope and slope variability in the land surface; it controls the probability of a contaminant running off or remaining on the surface in a specific area long enough for it to enter the aquifer 23 . The topography is described in terms of slope percentage 41 , and the topography slope map was constructed from the digital elevation model (DEM) using ArcMap10.4.1, with a cell size of 50 m × 50 m. The slope was then categorized according to the slope rating values, as shown in Table 7. Figure 5b illustrates the slope % for the study area.

Impact of vadose zone (I).
The vadose zone is defined as the zone above the water table that is unsaturated; it lies below the soil horizon and above the water table 30 . In the vadose zone, attenuation processes such as biodegradation, neutralization, filtration, chemical reaction, volatilization, and dispersion are likely to decrease as depth increases 23 . The rating for the vadose zone map was created using the groundwater flow system (GWFS) and the  www.nature.com/scientificreports/ material conformation for each of the aquifer types; conformation material is described in the Victorian aquifer framework 49 . Table 8 shows the vadose zone ranges and rating values proposed for the Pesticide DRASTICL model. Figure 5d shows the impact of the vadose zone map categorized with the rating values of the modified Pesticide DRASTICL model.

Hydraulic conductivity of the aquifer (C).
Hydraulic conductivity (C) is the aquifer's material capability to transmit water; it controls the flow rate of groundwater at a given hydraulic gradient, intergranular porosity, tectonic lineaments, and bedding planes 23 . Similarly, contamination is controlled by the flow rate of groundwater 41 . Furthermore, C controls the contamination movement in the aquifer 23 . Within index-based methods, C is denoted as the factor with the highest error associated with its estimation 60 . The factor rating for C was obtained from two sources. The first is an inferred value from various reports included in the GH GWFS by Dahlhaus et al. 48 The second is the data containing conductivity values from a groundwater transient model by SKM 45 . In this study, k values (m/d) were used from SKM 45 . A new rating scale is assigned to each range of hydraulic conductivity, as shown in Table 9. A rating map for C is shown in Fig. 5e.
Land use (L). Changes in land use induce the application of agricultural chemicals, industrial waste spills and landfill leachate, which can potentially leach into groundwater. L is a decisive and inducing factor of aquifer contamination thorough anthropogenic activities 60 . Research shows that specific methods that integrate information about land use perform better than intrinsic approaches 61 . Specifically, land use provides additional information on the use of the actual landscape, the type of cropping and crop rotation, which can be used to identify the amount and type of agricultural products used in the catchment 61 . This approach has been successfully applied in India 43 51 , and we selected the secondary land use cell of the Victorian land use information system (VLUIS) to assign new rating values from 1 to 9, as shown in Table 10. Figure 5f shows the land use rating map. Additionally, it shows the land use rating values for different land uses. Next, a raster file was created including the new ratings.
Multivariate statistical assembly of seasonal groundwater vulnerability. In standard DRASTIC methods, GWV has been estimated as a fixed or static value. However, such methods have failed to depict the impacts of temporal variations. Vulnerability can be affected by several factors, such as groundwater extraction, temporal variations in precipitation and evapotranspiration rates, temporal variations in groundwater depths, and temporal variations in groundwater recharge.
To make the vulnerability index sensitive to those changes, the proposed multivariate statistical analysis integrates time variations in depth and recharge while estimating seasonal GWV. This was achieved by estimating a GWV value for each season (spring, summer, autumn and winter). (2) I DSpr = Dr Spr * Dw Spr + Rr Spr * Rw Spr + Ar * Aw + Sr * Sw + Tr * Tw + Ir * Iw + Cr * Cw + Lr * Lw (3) I DSum = (Dr Sum * Dw Sum + Rr Sum * Rw Sum )+ Ar * Aw + Sr * Sw + Tr * Tw + Ir * Iw + Cr * Cw + Lr * Lw (4) I Daut = (Dr Aut * Dw Aut + Rr Aut * Rw Aut ) + Ar * Aw + Sr * Sw + Tr * Tw + Ir * Iw + Cr * Cw + Lr * Lw  www.nature.com/scientificreports/ After calculating the seasonal vulnerability index Eqs. (2)(3)(4)(5), we addressed the estimation of a new set of factors weights with the purpose for eliminating the correlation between factors 25 . This was achieved by separating DRASTICL factors into two hydrogeological groups: the dynamic group (D and R) and the static group (A. S, T, I, C and L). For both groups, an independent CA was performed, and a new set of upscaled factor weights was obtained for each season. CA deals with the subjectivity bias from using expert opinion for allocating factor weights.
This approach to data treatment (static and dynamic factors) is different from the standard index-based methods, which integrate all factor weightings in one multivariate analysis. CA has been selected because it has been shown to be the most adequate method for transforming the interrelated DRASTIC factors into uncorrelated DRASTICL vectors 25 .
First, the dynamic factors D, R plus a dummy factor (Q) with a constant value of 5 are resampled in ArcMap 10.4.1 using a fishnet at regular intervals of 5 km. The purpose of the sampling is to extract factor rating values at each point. Next, another point resampling is performed for each of the static factors A, S, T, I, C, and L. www.nature.com/scientificreports/    www.nature.com/scientificreports/ Second, a CA is performed to calculate vector loadings for each group. In the case of the dynamic group, a dummy dataset Q with a mean of 5 in rating value is used to assist in the calculation of the vector loadings as CA can only be performed with more than two factors. CA is also performed for the static hydrogeological group using the A, S, T, I, C and L factors. Estimating the time-variant Pesticide DRASTICL vulnerability. To improve the outputs of the Pesticide DRAS-TICL method, factor ratings were reduced to a maximum value of 9, and an independent multivariate CA was applied to both factor groups. After rating categorization, a resample is made using a fishnet of 5 × 5 km to obtain point values from each factor. The fishnet is the base for digital sampling of each of the factors included in the CA. The harmonization of scales method was applied to the resulting CA vector loading values, re-coordinating them from 1 to 5. As a result, seasonal Pesticide DRASTICL weighting values were used to estimate seasonal GWV.
The seasonal estimation of the GWV was quantified by using the corresponding seasonal factor weights. Eight raster layers from DRASTICL factors were used to calculate the seasonal I D . The results were mapped in the study area for each of the seasons.
Different values for vulnerability were expected for each season. The vulnerability values could be similar in places where the levels of depth and water recharge had small or no changes, indicating that the aquifer hydrogeological properties are less sensitive to climatic changes. Conversely, vulnerability values should be higher where depth and recharge showed large seasonal changes in depth and recharge.
Consent to participate. By this means, we the authors give explicit consent to participate in the submission of this scientific paper.

Correspondence analysis and time-variant groundwater vulnerability.
Vector loadings for the static groups L, A, S, T, I, and C are plotted in Fig. 6a, where the largest absolute value in loadings corresponds to S. Figure 6b shows the results for common vector loadings from the dynamic vector 1 D, vector 2 R and vector 3 Dummy (Dum) mapped for spring, summer, autumn and winter. www.nature.com/scientificreports/ The application of CA to data for the seasonal group resulted in the confirmation that the use of a dummy dataset (Dum) does not impact the value of the cumulative percentage of the system variance. It can be seen that Factor R has a larger loading compared to D. As shown in Fig. 6a, there are differences in the CA loading values for both groups; the larger the absolute value, the higher the non-correlation between factors. Such variables are called explicative factors; for the non-explicative factors, the values were the lowest 42 . In the case of the dynamic factor group, Factor R has a higher weight than Factor D. This could be explained because recharge is the driving hydrogeological process through which aquifer pollution occurs; the higher the recharge, the higher the chances for a pollutant to be transported to the water table 57 . Figure 6a shows the distribution of eigenvalues and the cumulative percentage of the system variance for common vectors V1, V2, and V3. Two explicative factors for the dynamic group vector 1 D and vector 2 R can explain 100% of the system variance. Figure 7b shows the loading values for V1 and V2 plotted for each season, and a clear difference in the loadings for R can be seen. Table 11 shows the cumulative percentage of the system variance and eigenvalues for each of the seasons. It can be seen that Dum has synthetic eigenvalues of zero and that both D and R can explain 100% of the system variance. Table 12 shows the results for dynamic CA factor loading values, which show that in absolute values, Factor R has higher loadings than Factor D.
The second CA is performed for the static group (Fig. 7a). Table 13 shows the results of the CA static eigenvalues and cumulative percentage of variance. Figure 8 shows the results of CA for the static factors, and common vectors V1, V2, V3 and V4 account for the cumulative percentage variance at approximately 98.6%. Figure 9 shows the distribution of static CA V1 and V2 loadings, which shows that vectors S, L and C are the explicative factors of the static system. Vectors S, L and C are the ones with the highest loadings and the highest distances from the centroid. Although the results from static CA suggest that principal static vectors V1, V2 and V3 can explain up to 93.7% of the system variance (see Fig. 8) and that using only these three factors should be enough to derive the vulnerability index, all static factors were considered for this study under the premise that each one is a significant explicative factor that contributes different types of information to the time-variant Pesticide DRASTICL vulnerability index.
New seasonal vector DRASTICL weights are presented in Tables 14 and 15 as a result of the harmonization of scales using the loading values for each of the factors. Such weights will be used in the calculation for the seasonal Pesticide DRASTICL vulnerability. The new vulnerability is referred to as the seasonal vector DRASTICL index and is represented as I DSpr , I DSum , I Daut , and I DWin .
From Fig. 10, it can be observed that there is a dimensionless scale for each season that should be approached as an independent representation of vulnerability. The results also showed that the risk was highly impacted by factors R, D, and L. The model was capable of estimating seasonal changes in vulnerability, with summer presenting the highest risk of contamination to groundwater in this specific area.   www.nature.com/scientificreports/   ) as NOX and dissolved reactive phosphorus (DSR) as phosphorus. Using the vulnerability map for winter, a field survey was designed and deployed at the study site for 2017 and 2018. Using the results from the model, 18 bores were sampled and analysed across the study site, as shown in Fig. 2. Figures 11 and 12 show the correlation coefficients (R 2 = 0.568 and 0.7056) for nitrogen and phosphorous, respectively. It can be observed that phosphorus outperforms nitrogen in evaluating the performance of the model. However, it should be considered that phosphorus showed fewer occurrences in groundwater and that the concentrations are lower compared to nitrogen. Researchers commonly use nitrogen to evaluate standard DRASTIC model performance.

Discussion
This study presents a novel approach to estimate a time-variant GWV as an alternative to a previously published static vulnerability 31 . Each seasonal vulnerability estimation is the representation of different conditions in the aquifer that reflects physical changes in the environment and should be interpreted accordingly. The resulting seasonal GWV is an independent estimation for each of the seasons. The qualitative categorization of vulnerability (very high, high, medium, low and very low), as mentioned in most previous research, may lead to the loss of valuable information for the user. In this study, a scale provides better insights into seasonal changes. The higher the vulnerability value on this scale, the higher the risk of contamination to groundwater.
Vulnerability values ranged from 42.25-179.89, 33.93-159.81, 34.08-168.74, and 45.56-205.20 for spring, summer, autumn, and winter, respectively. A lower vulnerability value was observed in summer, which is the season where less water is available for the aquifer due to less precipitation and high temperatures, producing the highest evapotranspiration rates and having a direct impact on the groundwater net recharge. Conversely, the highest vulnerability value occurred in winter, where the environmental conditions are the opposite of those in summer. This means that the vulnerability scale behaviour is consistent with the behaviour of the transient The seasonal CA suggests that the highest weight corresponds to the dynamic Factor R, followed by the static Factor L, while D showed a significant variation in weighting values between seasons, being highest in winter and spring. That is, in the season with the highest precipitation, D is weighted higher than in summer with the lowest weighting value.  www.nature.com/scientificreports/ In groundwater pollution surveying, it is common for researchers to use surrogate pollutants to predict the likelihood of pesticide pollution, as pesticide analysis is expensive and requires sophisticated laboratory facilities. Instead of analysing a large number of groundwater samples for pesticides, we selected a smaller number of samples that could show high concentrations of NOX and/or DRP for potential pesticide contamination. We also changed the original weights recommended by Allert et al. 30 using a correspondence analysis that derived new weights based on how sound a factor is compared to the others.
Phosphors showed a better correlation with the predictions of the model. A fundamental step in standard index-based methods is the calibration process 70 ; however, in this study, there are no seasonal data available for the seasonal calibration step, and this has been identified as a major gap requiring consistent effort and investment. We suggest that the calibration of the proposed time-variant Pesticide DRASTICL model should be carried out with seasonal nitrogen and/or phosphorus concentration data and that the continuous use and calibration of the model will lead to better performance.

Conclusions
In current index-based methods, the vulnerability index is mapped as a fixed scale assuming that such values do not change in the study area throughout the year. These methods ignore time variations in vulnerability. This study has addressed this knowledge gap by categorizing the influencing factors into two hydrogeological groups (static and dynamic) to estimate a time-variant groundwater vulnerability. In this approach, a Pesticide DRASTICL model was coupled with a double correspondence analysis to develop a seasonal vulnerability index.
The results of the model indicated that the vulnerability index varies between seasons and highly depends on the two dynamic factors, groundwater recharge and groundwater depth. We found that groundwater recharge and soil media are the explicative factors in the seasonal quantification of GWV. The model results showed a good correlation with the observed nitrogen and phosphorus levels. However, the model needs to be a recursive system of use, testing, and calibration, as there is not enough seasonal data. The quantification of temporal and spatial variations in the GWV will provide additional information to farmers and water administrators to better manage groundwater pollution prevention programs. It will also assist with the use of fertilizers and pesticides in agricultural areas.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The code that supports the findings of this study is available from the corresponding author upon reasonable request.